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Abstract 

The supersymmetric evolution of the fragmentation functions (or timelike 
evolution) within N = 1 QCD is discussed and predictions for the fragmenta- 
tion functions of the theory (into final protons) are given. We use a backward 
running of the supersymmetric DGLAP equations, using a method developed 
in previous works. We start from the usual QCD parameterizations at low 
energy and run the DGLAP back, up to an intermediate scale -assumed to 
be supersymmetric- where we switch-on super symmetry. From there on we 
assume the applicability of an N = 1 supersymmetric evolution (ESAP). We 
elaborate on possible application of these results to High Energy Cosmic Rays 
near the GZK cutoff. 
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1 Introduction 



The Standard Model of elementary particle physics successfully accounts for all the 
observed particle properties. Its remarkable success suggests the possibility that 
it correctly describes the particle properties up to energy scales which far exceed 
the currently accessible energy scales. Indeed, this possibility has been entertained 
in the explorations of Grand Unified, and Superstring, Theories. Additionally, the 
spectacular confirmation of the Standard Model in high energy colliders, as well as 
the proton longevity and suppression of left-handed neutrino masses, provide strong 
support for the grand desert scenario and unification. The question then is how are we 
going to find out whether this is indeed the path chosen by nature, as direct signatures 
at LHC/NLC will be able to probe the desert only up to a few TeV. Therefore we 
need new tools to achieve this. An example of such a probe that has been used with 
great success is that of the gauge coupling unification. The general methodology in 
this respect is to extract the gauge couplings from experiments that are performed at 
the accessible energy scales. Based on concrete assumptions in regard to the particle 
content in the desert, the couplings are then extrapolated to the high energy scale and 
the unification hypothesis is tested. In this manner the consistency of the specific 
assumptions in regard to the particle content in the desert with the hypothesis of 
the gauge coupling unification is subjected to an experimental test. The well known 
spectacular success of this methodology is in differentiating between gauge coupling 
unification in supersymmetric and non-supersymmetric Grand Unified Theories [|T|. 
The unification hypothesis is found not to be consistent with the low energy data, 
unless the minimal particle content is supersymmetric, or is modified in some other 
way. Our task then is to develop additional tools that can serve as similar useful 
probes of the desert. 

1.1 Motivations for this work 

In this paper we start applying this philosophy to another physical setting. 

A different experimental observation - which also points toward the validity of the 
big desert scenario - is the existence of Ultra High Energy Cosmic Rays above the 
Greisen-Zatsepin-Kuzmin cutoff. A plausible explanation for the observation of such 
cosmic rays is the existence of supermassive metastable states with mass of the order 
10 12 — 10 13 GeV, and a lifetime that exceeds 10 w y. The possibility that such states 
compose a substantial component of the dark matter as well as that they may explain 
the observation of UHECR has been discussed elsewhere. It has been further shown 
that realistic string models often give rise to states with precisely such properties 0. 

The existence of such super-massive, slowly decaying states, may therefore serve 
as a probe of the grand desert. A plausible explanation of the observed showers in 
the UHECR is that the original supermassive decaying particle decays into strongly 
interacting particles which subsequently fragment into hadrons ||. This hypothesis 
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therefore needs the relevant fragmentation functions at the energy scale of the origi- 
nally decaying particle. These functions, however, similarly to the case of the gauge 
couplings considered in much of the literature on unification, are extracted from ex- 
periments at the currently accessible (low) energy scales which are much below that 
of the decaying particle. It is obvious that this points toward the use of the renormal- 
ization group (RG) evolution. However, this evolution, as well as the composition of 
the produced shower, depends on the assumptions made in regard to the composition 
of the spectrum in the desert. We conclude that similarly to the gauge coupling uni- 
fication methodology, the evolution of the fragmentation functions by utilizing the 
renormalization group extrapolation, and the subsequent implementation in the anal- 
ysis of the UHECR hadronic showers, can serve as a tool that differentiates between 
different assumptions on the very high energy spectrum. 

With these motivations in mind, here we develop the instruments to evolve the 
fragmentation functions in supersymmetric QCD up to a high energy scale. In par- 
ticular, in this paper we present for the first time the structure of the supersymmetric 
DGLAP (Exact Supersymmetric DGLAP or ESAP) equations in the timelike region. 
We solve them numerically for all flavours and introduce various combinations of non- 
singlet matrix equations to reach this result. We then study the impact of the new 
evolution assuming a common scale for susy breaking and restoration, parameterized 
by the mass of the superpartners, here assumed to be mass degenerate. 

1.2 Step Approximations 

Given the long stretch it takes to proceed with the analysis of the susy evolution, a 
topic of numerical complexity in its own, the applications to cosmic rays of our work 
will be presented in a companion paper that will follow shortly. We also have decided 
to address issues related to the fragmentation region of N — 1 QCD starting from 
the "low energy" end, since nothing is known of these functions at any scale. The 
term "low energy" (the results we discuss here are obtained in the 10 3 — 10 8 GeV 
energy range ) - from the point of view of collider phenomenology - is a misnomer, 
but it is not in the context of high energy cosmic ray physics. There is also another 
more direct reason for limiting our analysis to this range. We have discovered some 
features of the RG evolution which require a special care at such large energies. We 
find -in our numerical studies- the onset of an instability in the evolution in the 
singlet sector of the fragmentation which requires an independent investigation. It 
is not clear to us, at this point, whether this instability is an intrinsic limitation of 
the algorithm used in the numerical implementation or if the perturbative expansion 
needs a resummation. At such large scales this last option remains open. We should 
also mention that we work in "x" space, and this might be a limitation. We hope to 
return to this issue in the near future. 

The analysis of fragmentation functions of Susy QCD is a topic of remarkable 
phenomenological interest both for collider phenomenology and in the astroparticle 
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context. While - recently - a detailed analysis of the supersymmetric evolution of 
the parton distributions of Susy QCD has been presented || |[ ||, the study of the 
evolution of the supersymmetric fragmentation functions is still missing. Of partic- 
ular interest is the study of a combined QCD/SQCD evolution with intermediate 
regions in the evolution characterized by partial supersymmetry (SAP) or exact su- 
persymmetry (ESAP) ||. It is well known that a matching between these regions is 
possible with specific boundary conditions. These boundary conditions should come 
from some deeper understanding of the way in which supersymmetry is broken and 
restored as we run into the different stages of the renormalization group evolution. 
In a first approximation, however, it is possible to assume that the regular QCD dis- 
tributions (or fragmentation) functions are continuous at each intermediate region, 
thereby neglecting threshold effects which can't be inferred from first principles. As 
we cross any supersymmetric region, starting, for instance, from the end of the reg- 
ular QCD evolution, it is possible to assume that supersymmetric distributions and 
fragmentation functions are generated radiatively. We recall that a similar modeling 
-which neglects possible effects of higher twist at the thresholds- is well spread in the 
case of ordinary QCD. 

2 Initial/Final State Scaling Violations 
2.1 The Spacelike Evolution 

The first analysis of the spacelike evolution was carried out in ref. [[?]], where simple 
models -obtained from the analysis of the first 2 moments - of the supersymmetric 
parton distributions were also presented. A complete strategy to solve these equations 
and the strategy to generate supersymmetric scaling violations in N=l QCD was 
put forward in Refs. j|, In N = 1 QCD gluons have partners called gluinos 
(here denoted by A) and left- and right-handed quarks have complex scalar partners 
(squarks) which we denote as q~L and q~ R with q = qi + qn (for left-handed and 
right-handed squarks respectively). 

The interaction between the elementary fields are described by the 577(3) color 
gauge invariant and supersymmetric Lagrangian 

+iqi pq { + D^q R D»q R + D^q L D»q L + igsft (\ a R q\ L T a q Li + \ a L q\ R T a q Rl - h. c.) 
-\d 2 (q[i Ta QLi - q Ri T a q Ri ) 2 + mass terms, 

(1) 

where a runs over the adjoint of the color group and % denotes the flavours. 
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Figure 1: An illustration of the regular versus a mixed supersymmetric evolution 

qiR are the supersymmetric partners of right-handed quarks qm and qn are those 
of the left-handed quarks qn. 

We remind here that in actual phenomenological applications, it is useful to ex- 
tract the light-cone dynamics -or parton model picture- of supersymmetric collisions 
by invoking a factorization of the cross section into hard and soft contributions. Frag- 
mentation functions appear in this context. An integral part of this analysis is the 
use of renormalization group equations -susy DGLAP- which resum the logarithmic 
scaling violations and evolve distributions functions and fragmentation functions to 
the appropriate factorization scale. 

These studies require the matrix of the anomalous dimensions (for all the mo- 
ments), or of the corresponding Altarelli-Parisi (DGLAP) kernels. We recall that the 
leading order anomalous dimensions are known |8], |7j both for a partial susy evolution 
and for the exact one. 

2.2 The Timelike Evolution 

In this work we focus our attention on the evolution of fragmentation functions in the 
context of exact supersymmetry, with coupled gluons and squarks. Fig. 1 summarizes 
the basic strategy of our work. The figure at the top depicts a regular evolution (a 
one-phase evolution), while the one at the bottom illustrates a mixed QCD/SQCD 
evolution (a two-phase evolution). The first stage is supersymmetric, the second 
one is regular. In the first stage all the partons (susy and not susy) are massless, 
in the second one only quark and gluons survive. At the end of the first stage the 
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super symmetric partners are extinct, but their presence at the higher scales is saved 
in the boundary conditions for the quarks and gluons fragmentation functions, when 
the second stage of the evolution starts. 

Within the approximation implied by the use of a pure SQCD evolution, no 
missing E t events are allowed along the development of the supersymmetric cascade, 
since electroweak and flavour mixing effects are not included. Their inclusion is still 
an open chapter even in (QCD) Monte Carlo event generators for the final state. 

The ladder structure of the supersymmetric showers in the leading logarithmic 
approximation can be easily pictured in terms of symmetric doubling of lines of 
regular QCD ladders in all the possible allowed cases. Showers can be equally initiated 
by q q pairs, gluon pairs, gluino pairs or squark pairs. The supersymmetric transition 
(i.e. those involving supersymmetric partners) stop (see Fig. 1 (bottom)) once the 
supersymmetry breaking scale M\ = Mq is reached. 

Similarly to the QCD case, in the case of exact N — 1 supersymmetry we de- 
fine singlet and non-singlet fragmentation functions Dj(x,Q 2 ). They describe the 
amplitude for a parton of type / to fragment into a hadron h as a function of the 
Bjorken variable x (fractional energy of the fragment) and initial energy Q. Their 
operatorial definition is similar to that of ordinary parton distributions. The evolu- 
tion equations for these functions are similar to the standard DGLAP equations, but 
with transposed kernel matrix (P — ► P T ) . In leading order the analytic continuation 
(from the spacelike or DIS evolution to the timelike evolution) is straightforward, 
while some complications appear to next-to-leading order. They involve a breaking 
of the Drell-Yan-Levy relation, which appears to be violated at parton level |J. Since 
our analysis is in the leading logarithmic approximation, we will not consider these 
aspects in our work any further. 

The equations for the timelike evolution are given by 
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rflog(Q 2 ) 



+Pm®Dl) (5) 



(6) 



In the following we will use the short-hand notation 

Q 2 ) = D^(x,g 2 ) + Dj 4il (x,g 2 ) (7) 

to denote the fragmentation functions of squarks of flavour % at a fractional energy 
x and momentum Q. It is also convenient to separate the equations, as usual, into 
singlet and non singlet sectors using the definitions 

n f 

D h qv (x,Q 2 ) = EK(^ 2 )-^,(^ 2 ))- 
i=i 

= D h 

QNS 
nj 



dim 2 ) = EK(^V4,(^ 2 )) 

1=1 

- D LM 2 ) 

D q+ (x,Q 2 ) = £(D h qi (x,Q 2 ) + D h qi (x,Q 2 )) 
i=i 
n f 

/jj + (x,g 2 ) = ^(/j|(x,g 2 ) + /j|(x,g 2 )). (8) 
i=i 

The non singlet equations are 



dQ 2 ' 2tt 

Q 2 ^ 2 D qv M 2 ) = ^(P q «®D h qv + P m ®D^), 



(9) 



and the singlet matrix equations, which mix qy and q v with the gluons and the 
gluinos 
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D\(x,Q 2 ) 

Dq+{X,Q 2 ) 
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P 

1 gq 


P - 1 
1 gq 




Pxx 


Pxq 


Pxq 


P 

r qg 


Pq\ 


P 

r qq 


P - 

r qq 


p. 
L r qg 


PqX 


P- 
r qq 


P~~ 



T 



D h Jx,Q 2 ) 
D\{x,Q 2 ) 

D h 5(+) {x,Q 2 ) 



(10) 



where "T" indicates the matrix transposed. To solve for all the flavours, it is conve- 
nient to introduce the linear combinations 



Xl (x,Q 2 ) = D\ +) 



n f 
1 



D h qW 



x,(x,g 2 ) = D h (+) -— D h - M 
% n f 



(11) 



and the additional singlet equations 



(-) 



dQ 



2 5, 



(-) 



%Q 2 ) = 
%Q 2 ) = 



a(Q 2 ) 

a(Q 2 ) 
2tt 



qq 



qq 



^ V) + Pqq 



^(-) + Pqq 



D h (-) 



(12) 



and 



2, ( P - 



2tt 
a(Q 2 



qq 



qq 



D + P- 
u x l + r qq 



D i + Pi 



qq 



(13) 



The general flavour decomposition is obtained by solving the singlet equations for 
Dy (+) and D~ M , then solving the non-singlet equations for D\_) and D\ +) and for 

D } 1 and D% . The fragmentation functions of the various flavours are extracted using 
the linear combinations 



Dt 

Hi 



Hi 



1 

2 



D\-) + D h Yi + -D% 



■ ' n s ~ qy 



- 1 - U. 

2 V q * 



1 

n f 



Dl - -D h q(+) 



(14) 
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and 

(15) 

to identify the various flavour components. 

There are simple ways to calculate the kernel of the susy DGLAP evolution by 
a simple extension of the usual methods. The changes are primarily due to color 
factors. There are also some basic supersymmetric relations which have to be sat- 
isfied ||. They are generally broken in the case of decoupling. We recall that the 
supersymmetric version of the f3 function is given at two-loop level by 



f3% = -(llCU-2n/-27i A ) 

P? = ^(3AC A -10C A n f -10C A n x -6C F n f -6C x n x ) (16) 

where n/ is the number of flavours and for Majorana gluinos n\ = Ca- 

The ordinary running of the coupling is replaced by its supersymmetric running 

* S (Ql) = 2 1 / _ /3f lnln(Q 2 /A 2 ) 1 A 

2vr ~ /? 5 ln(Q 2 /A 2 ) ^ /f ln(Q 2 /A 2 ) + l ln 2 (Q 2 /A 2 )V 1 ) 

The kernels are modified both in their coupling (a — > a s ) and in their internal 
structure (Casimirs, color factors, etc.) when moving from the QCD case to the 
SQCD case. In order to illustrate the approach that we are going to follow in our 
analysis of the fragmentation functions within a mixed SQCD/QCD evolution, we 
recall the strategy employed in Refs. j|, [5] to generate the ordinary supersymmetric 
distribution functions. 

In Ref. |J for scaling violations affecting the initial state were introduced 3 regions: 

• 1) the QCD region, described by ordinary QCD (Altarelli Parisi or AP) ; 

• 2) an intermediate supersymmetric region with coupled gluinos and decoupled 
squarks (partially supersymmetric AP or SAP); 

• 3) the N=l region (exact supersymmetric AP or ESAP). 

The basic idea was to generate these distributions radiatively as usually done in 
QCD for the gluons, for instance, using the fact that the matrix of the anomalous 
dimensions is not-diagonal. The most general sequence of evolutions (denoted AP- 
SAP-ESAP in |Q) is described by the arrays {Qi,Qf)Ap — {Qi,Qf)sAp — {Qi,Qf)ESAP, 
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with Qf,AP = Qi,sap = and Qj,sap = Q%,esap = iriq- In this work we limit our 
analysis to a simpler AP-ESAP evolution and we run the DGLAP equations back 
starting from known QCD fragmentation functions - for which various sets are avail- 
able in the literature fl ij - up to an intermediate supersymmetric scale. From this 
point on (up in energy) we switch on the ESAP evolution of fragmentation functions. 
As we run the equations upward, supersymmetric fragmentation functions are gener- 
ated. The boundary values of the low energy (QCD) functions at the supersymmetry 
scale m 2 A set the initial condition for the (backward) supersymmetric running up to 
the final scale. In simple terms: we reach the mountain from the valley, and cross a 
fence along the way. 

In general, if we split the intermediate susy scale into 2 sectors m q << m^x 
(ESAP-SAP-AP evolution) in this (general) case the solution is built by sewing the 
three regions as 



Dl(x,Q 



D h q Jx,Q 2 ) 



+ I dlog Q 2 Pl> P NS (x,a(Q 2 )) 

JQ.i 



DL(x,Q 2 



+ 



+ 



dlogQ 2 P^ A N P s (x,a^Q 2 ))® 



f Qf d\ogQ 2 P^ s N A s P (x,a ES (Q 2 ))® 

J mi. 



D h qv {x,Q 2 ) 


D h qv (*,Q 2 ) 
Dl(x,Q 2 ) 



in the non singlet and 



;is) 



D h G (x,Q}) 
D h x (x,Q 2 ) 
D h + {x,Q 2 ) 



+ 



+ 



D h G (x,Q 2 ] 



x, Q 2 e 



2 

log Q 2 P» s AP (x,a(Q 2 ))® 



d\ogQ*P$ s SAP {x,a b {Q''))® 



x,Q 2 ) 



D h G ( 

D h x (x,Q 2 ' 



D 



(x,Q 2 ) 




\ Q ] d\ogQ 2 P^I SAP (x,a ES (Q 2 ))® 



D h G (x,Q 2 ) 
D h x (x,Q 2 ) 
D h q+ {x,Q 2 ) 
D~+(x,Q 2 ) 



D h G (x,Q 2 



D h 





+ \x,Q 2 ) 




(19) 



for the singlet solution. The zero entries in the arrays for some of the distributions are 
due to the boundary conditions, since all the supersymmetric partners are generated, 
in this model, by the evolution. 
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The general structure of the algorithms that solves these equations is summarized 
below. We start from the non singlet sector and then proceed to the singlet. From now 
on, to simplify our notation, we will omit at times the index u h" from the fragmen- 
tation functions when obvious. We define D q Ns(x,Q 2 ) = (D qv (x,Q 2 ), D qv (x,Q 2 )) T 

and DA n (x) = (jDA^ , and introduce the ansatz 



D qNS (x, Q 2 ) 



n=0 



71! 



(20) 



where n is an integer at which we stop the iteration, which usually ranges up to 20. 
The first coefficient of the recursion is determined by the initial condition 



DA (x) = U(x)®D" s (x,Q 2 ), 



(21) 



where 

Ui(x) = 

U 12 (X) EE 

The recursion relations are given by 

DA n+1 (x) = - 



5(1 -x) 



8(1 -x) 

5(1 -x) 



(22) 



Pi 



pT,NS 



DA n (x) 



(23) 



The solution in the first (DGLAP) region at the first matching scale m 2 \ is given by 

(24) 



no 



sr ^DA AP (x) 1 n (a(m 2X 2 )' 

C ,„ s( ,, m2A ) = E _^log»(-L^ 



n=0 



At the second stage the (partial) supersymmetric coefficients are given by (S is a 
short form of SAP) 



DA^ NS (x) = U(x)®D q (x,m 2X 2) 



DA^ s (x) = ~P^(^Mf s (4 
Po 



(25) 



We construct the boundary condition for the next stage of the evolution using the 
intermediate solution 



(26) 
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evaluated at the next threshold m 2 \- 



D ^ ) = t£mU^L). (27) 



n=o n - 



a(m 2 \ 2 ) ) ' 

The final solution is constructed using the recursion relations 

DA ES (x) = U(x)®q(x,rriq) 
DA E n s +1 {x) = -^r s P%l{x)®DA E n s {x) 

The final solution is written as 



(28) 



n ( rfi\ n < ^UV^i f ^m 2X 2 ) \ " DA^jx) ( a s (m}) ' 
D M x,Q) = ^(,,Qo) + E^^log^^^J + E^^log^^^ 

"° DA ES (x) ( a ES (Qj) \ 



+ e iQ g 



71=1 



n! la^^m^ 2 ) 



(29) 



As we have already mentioned above, in the analysis presented below the two susy 
scales will be collapsed into one (m 2 \). 

3 QCD Evolution of the Fragmentation Functions 

As we have already mentioned, the timelike and the spacelike evolution, in leading 
order, are essentially the same. The singlet kernels are just the transposed of the 
kernels describing the evolution of the parton distributions. The ordinary QCD 
kernels, in leading order, are given by 



P (o) _ r 



1 +X 



2 % 



qq,NS " °F h 



X 



p(o) (r \ _ p(0) 

r qq \ X ) — r qq,NS 

P£\x) = n, (x 2 + (1 - xf) 

pSHx) = c F l - + ( l - x? 



x 



P%\x) = 2N c (^-^ + ±-2 + x(l-x))+^5(l-x) (30) 
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with Pq = 11/3 Ca — ^/STnrif being the first coefficient of the QCD /3-function 
and Tr = 1/2. rif is the numbers of flavours. The equations for the fragmentation 
functions in QCD are given by 



d\og(Q 2 ) qiK ' 2tt V 2n/ 



2^ 



and, as usual, can be decomposed into a non-singlet and a singlet sector 



(31) 



g2 d^V )(x ' g2) = ?^P qq ,N S (x,a(Q 2 ))®D h qt) (x,Q 2 ) (32) 

for the non-singlet distributions and 

d ( D h ai+) (x,Q 2 ) \_( P qq P 9q \^( D%{x,Q*) \ 



H 










P 99 ) 



d\og(Q*)\ D*(x,Q 2 ) J \P qg P gg )^\ D*(x,Q*) 
for the singlet sector, where 



(33) 



^ ( _,(x,g a ) = D* -D*. (34) 



A fast strategy to solve these equations, as discussed in |4j], where a complete leading 
order evolution has been implemented, is to solve the recursion relations numerically. 
In the timelike case -that we are considering- these relations are obtained as linear 
combinations of the spacelike ones. We get 



,(-) 



n i". a Cf ( 1 dy yDA q T\y)-xDA c t\x) C F m , 5 H, , 

DAn+i = -4-5-/ 4— log(l -x)DA% (x) 

Po Jx y y-x Po 



Po J* y Po 



^ = _ 4 g£ /■ jg^-'M - !gcw _ 4 g£, og(1 _ , )JMrH 

Po Jx y y — x p 



+2 £L f*( 1+l|M f> W -3| M f> (l ) 
Po Jx y po 



(35) 
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A similar expansion is set up for the non-singlet variable Xi — ~ V n / 

DA% +1 = -4% f *L V DA ?M ~_ * x DA *M - £JL log(i - X) DA?( X ) 

+2 Cf j ld JL {l + z) DA x i{y) _ 2°L D A*t{x). 
Po Jx y po 

(36) 

Similarly, the singlet equations generate recursion relations of the form 

Po Jx y y — x po 

+ ^/;f(i + ,)^%)-3|^'(») 

_ 2 | /; ?ii±(±pL DA , (y) 

DAUM - -4 f - XDAl(X) - 4 W - *>^M 

Po Jx y y — x po 

_2 I f 7 + (1 " I)2 ) M: " ,,!,) - 

- 4 ir7(i- 2+;(i - 2) )« < 37 » 

The coefficients of the various flavours both for D„ and D% are then obtained 
from the relations 

£M* = ^IM* - ^4* (_) ) (38) 

4 The ESAP (N = 1 QCD) Fragmentation 

Moving to N — 1 QCD, we introduce recursion relations for appropriate linear com- 
binations of non-singlet fragmentation functions 



DAS* = rg ^M-^W -4^1og(l-,)lM» W 

% Jx y y-x Po 



Cp r^dy yDA% (y) - xDA* (x) „ C F 

A) ^ 1/ 

f^'M - 2^ /' ^zD^(„) + 2?f[ 1 ^(l + z ) DA*<y) 
Po Po y Po J* y 

(39) 
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-2^DA^{x) - 2^f f ^-zDA%(y) + 2^f f ^ (1 + z) DA%(y) 
Po Po Jx y p Jx y 

(40) 

and two similar non-singlet equations for D h ,_) and D h ( s 

-2^W-'(,) - 2^ /' *^-(v) + 4 /' * (1 + *) D^-'tv) 
Po Po A y Po Jai y 

(41) 



DA& = -4^ r ^ ggjtSO - ggjTw _ ^ lQg(1 _ x)d4%) 
p J x y y — x Po 

-2-fDA% (x)-2-f -^zDA% (y) + 2-f -H- (1 + z) DA% (y). 
Po Po Jx y Po Jx y 

(42) 

Analogous equations are satisfied by the combinations D^-) and D ! } ( _ ) by replac- 
ing in (H) and §|) L>V) D g(-) and £>-(-) D ^)- 
Moving to the singlet sector we obtain 



DAUM = _ 4 §1 /' ^ y£mH - ^M£> _ 4 £l log(l - X )DAI( X ) - DA i(x) 
Po Jx y y-x Po 

+ 2^d /' * (1 + ,) DA-(y) - f d JL (1 + , _ 2 ) M » W 
Po Ay po Jx y \z / 

+2 §i /■ ^ + (i - ,)•) !,>!.(,,) -2^ r^c^+d- 

Po Jx y ^ ' Po Jx y v 7 

_ 2 Ca ^ ^ + (1 _ ^ DA , (y) _^L( ld JL z{l _ z) DAUy) 
Po Jx y po Jx y 

(43) 

Po Jx y \z J Po Jx y y-x 
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-£± log(l - x)DA&x) + 2§i f ^ (1 + z) DAl{y) - DA*i 
Po Po Jx y 

Po Jx y Po 



x) 



_ 2 r}± f l djj (1 _ z) DAl(y) _ ^ ^ zDAi[y) (44) 
Po Jx y Po Jx y 



DAt + \{x) 



C F /-i dy /2 
Po •/* y \z 

(+) , 



2) DAi(y)-2^f ( ld f{l-z)DAl(y) 
/ po y 



Cf ^dyyDAtr'{y)-xDAi 
Po 



,(+), 



2/ 

(/.;/ 



7 (+) , 



V 

4^1og(l - x)DA( + \x) 
Po 



+2^ / ^ (1 + z) DAf (y) - 2^ZMf ) (x) - 2^ / ^zl^ (y) 
Po •/* y Po Po Jx y 

(45) 



dy 



IMfifo) = -2^ /'*(*- 2) DAi(y) - 2^ /' * (dA^(x) - DA^(y)) 
Po Jx y \z J Po Jx y ^ 1 

_ 4 CV ^dy yDAi(y)-xDAj(x) | f F f l dy^- 

Po Jx y y-x Po L y 

-4^r Ml " *)^f } (v) " ^DAf ) (x) (46) 
Po Po 

where z = x/y. 

5 Numerical Results 

As an illustration of the procedure we adopt in our studies, let's consider the decay of 
a hypothetical massive state of mass 1 TeV into supersymmetric partons. The decay 
can proceed, for instance, through a regular qq channel and a shower is developed 
starting from the quark pair. The N — 1 DGLAP equation describes in the leading 
logarithmic approximation the evolution of the shower which accompanies the pair, 
and we are interested in studying the impact of the supersymmetry breaking scale 
(m\) on the fragmentation. In our runs we have chosen the initial set of Ref. flTT 



The parameterizations are shown in the appendix just for illustrative purposes. In 
our analysis we focus on proton fragmentation functions. As in Ref. [11]], we introduce 
the scaling variable 

a = In ) ' ■ (47) 
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In LO we have Ag, = 88 MeV JTTJ. We use three different values for /x , namely |Tl| 



y/2 GeV, if a = u,d,s,g 

^ = <{ m (T) c ) = 2.9788 GeV, if a = c . (48) 

m(T) = 9.46037 GeV, if a = b 

This leads to three different definitions of s. For definiteness, we use the symbol s c 
for charm and for bottom along with s for the residual partons. We parameterize 
the f.f.'s as 

D(x,fi 2 ) = Nx a (l -xf (l + -) (49) 



x, 

and express the coefficients N, a, /3, and 7 as polynomials in s, s c , and s&. For 



s = s c = Sf, = 0, the parameterizations agree with Eq. (2) of Ref. JTT| in combina- 
tion with the appropriate entries in Table 2 of that paper. The charm and bottom 
parameterizations must be put to zero by hand for s c < and s& < 0, respectively. 

Typical fragmentation functions in QCD involve final states with p, p, 7r ± , ir° and 
kaons k ± . We have included the relevant sets from ref. |TT[ that we have employed 



in our numerical analysis in an appendix. We refer to the original literature for 
a list of all the fragmentation sets. We have chosen an initial evolution scale of 
10 GeV and varied both the mass of the susy partners (we assume for simplicity 
that these are all degenerate) and the final evolution scale. Since our concern is in 
establishing the impact of supersymmetric evolution and compare it to standard QCD 
evolution across large evolution intervals, we plot the initial fragmentation functions, 
the regularly evolved QCD functions and the SQCD/QCD evolved ones. The latter 
two are originated from the same low energy form of Ref. [|ll]] . The figures show that 
the effects of supersymmetric evolution are small within the range described by the 
factorization scales Qf and Qi (Qf = 10 3 GeV, Qi = 200 GeV). We mention that 
Qf is the starting scale (the highest scale) at which the decay of the supersymmetric 
partons starts. Qi is fixed by the gluino/squark masses and coincides with them. It 
is possible to include in the evolution of the fragmentation functions also different 
thresholds associated with more complex spectra in which the susy partners are not 
degenerate. These effects are negligible. Threshold enhancements may require a 
more accurate treatment and will be discussed elsewhere, however we don't expect 
them to play any important role, especially since we are interested in very extended 
renormalization group runnings. Fig. 2 shows the initial condition for the up quark 
fragmentation function. We have chosen Q = 10 GeV, extrapolated from collider 
data in Ref. JTT]] in their global analysis. The evolution of this function follows 
QCD from this lowest scale up to a scale of 400 GeV (th^a)- Above this scale we 
use the full N — 1 evolution. The initial fragmentation scale, here denoted by Qf 
is 1 TeV. In general, the small-x (diffractive) region gets slightly enhanced. The 
highest scale is not large enough to allow a discrimination of susy effects from the 
non supersymmetric ones. The trend of the two evolved curves is similar. In Fig. 3 an 
analogous behaviour is found for the fragmentation functions of charm and strange 
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quarks. Again, small-x enhancements are seen, but regular and susy evolution are 
hardly distinguishable. The situation appears to be completely different for the gluon 
fragmentation functions (f.f's) (Fig. 4). The regular and the SQCD evolved f.f.'s differ 
largely in the diffractive region, and this clearly will show up in the spectrum of the 
primary protons if the decaying state has a supersymmetric content. As we raise 
the final evolution scale we start seeing more pronounced differences between regular 
and supersymmetric distributions. We have shown in Fig. 5 the squark f.f.'s for all 
the flavours and the one of the gluino for comparison. The scalar charm distribution 
appear to grow slightly faster then the remaining scalar ones. The gluino f.f. is still 
the fastest growing at small-x values. 

In Figs. 6 and 7 this trend is quite evident for the quark and gluon fragmentation 
functions respectively. The latter is the distributions which remains by far the most 
sensitive to changes in the highest fragmentation scale. It is therefore very likely that 
the spectrum of final protons is substantially modified by this fast modifications of 
the gluon density. A sustained stability -as we raise the supersymmetric scale - of all 
the distributions is also noticed. In fact, in Fig. 8 we vary the gluino mass for a fixed 
final highest scale and notice a mild variation of the fragmentation function of the up 
quark. A similar behaviour is recorded in Fig. 9, where we compare the fragmentation 
functions of the up quark and of its superpartner at various m A 's. The shape of the 
gluino density as a function of the the initial fragmentation scale is shown in Fig. 10. 
The growth appears to be fast and is rather pronounced. Fig. 11 finally illustrates 
the behaviour of f.f's of all types at a large initial fragmentation scale. In order 
to illustrate the different behaviours that the quark/squark sector has compared to 
the gluon/gluino sector, we have shown the f.f.'s for the b quark and the b squark. 
Although all the functions get a small-x enhancement or are more supported in this 
x-region, the largest fragmentation appears to come from the gluon-gluino sector. 

6 Conclusions 

We have presented a first detailed study of the fragmentation functions of susy QCD 
in the leading logarithmic approximation and presented, for the first time, results for 
the f.f.'s of all the partons into protons within a radiatively generated model. Our 
analysis is motivated by the growing interest that high energy cosmic rays experiments 
will be receiving in the next several years from the astroparticle/high energy physics 
community. There are hints, from this work, that supersymmetric effects may affect 
substantially the fragmentation region and the spectra of the hadronic component 
of cosmic rays which is expected to play a dominant role as one approaches the 
GZK cutoff. Our analysis indicates that fragmentation probabilities are redistributed 
among the various open channel and have an impact on the spectrum of the primaries. 
Compared to initial state scaling violations where only a gluino of small mass seems 
to enhance the radiative structure of parton distributions [|J ^j, the cosmic ray 
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story clearly points to another direction. Given the very large energy available in 
the fragmentation region, and assuming a supersymmetry breaking scale of the order 
of 1 TeV, the f.f 's, which control the multiplicities of the hadron component of the 
various channels are substantially affected by the presence of a supersymmetric scale 
at the highest energy. The natural question to ask is then: what are the distributions 
of supersymmetric partons in a decaying metastable state of very large mass prior 
to fragmentation? At such large scales, so far from the 1 TeV scale preferred by 
many standard supersymmetric arguments, the Renormalization Group strategy can 
be easily embraced to its fullest extent with sizeable consequences on the low energy 
end. On this and other related issues we hope to return in more detail in the near 
future. 
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7 Appendix 1. Input Functions 



For convenience we have included below the list of f.f.'s taken from Ref. |TT[ that we 
have used 

N = 0.40211 - 0.21633s- 0.07045s 2 + 0.07831s 3 
a = -0.85973 + 0.13987s- 0.82412s 2 + 0.43114s 3 
P = 2.80160 + 0.78923s- 0.05344s 2 + 0.01460s 3 
7 = 0.05198s - 0.04623s 2 (50) 



Dp/pi 



2\ 



N = 4.07885 - 2.97392s- 0.92973s 2 + 1.23517s 3 

a = -0.09735 + 0.25834s- 1.52246s 2 + 0.77060s 3 

P = 4.99191 + 1.14379s- 0.85320s 2 + 0.45607s 3 

7 = 0.07174s - 0.08321s 2 (51) 



DP/P(x^ 2 ): 



N = 0.11061 - 0.07726s c + 0.05422s 2 - 0.03364s: 
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a = -1.54340 - 0.20804s c + 0.29038s 2 - 0.23662s 3 

P = 2.20681 + 0.62274s c + 0.29713s 2 - 0.21861s 3 

7 = 0.00831s c + 0.00065s 2 (52) 



N = 40.0971 - 123.531s 6 + 128.666^- 29.1808^ 

a = 0.74249- 1.29639s 6 - 3.65003s 2 + 3.05340s 3 . 

P = 12.3729- 1.04932s 6 + 0.34662s 2 . - 1.34412s 3 . 

7 = -0.04290s 6 - 0.30359s 2 (53) 



N = 0.73953 - 1.64519s + 1.01189s 2 -0.10175s 3 

a = -0.76986 - 3.58787s + 13.8025s 2 - 13.8902s 3 

P = 7.69079 -2.84470s - 0.36719s 2 - 2.21825s 3 

7 = 1.26515s- 1.96117s 2 (54) 

8 Appendix 2. The Weights of the N=l Kernels 

We briefly recall the numerical strategy employed in this analysis. A more detailed 
description will be given elsewhere. We just mention that the radiative generation of 
super symmetric distributions requires special accuracy since these scaling violations 
grow up very slowly. We define P(x) = xP(x) and A(x) = xA(x). We also define 
the convolution product 

The integration interval in y at any fixed x-value is partitioned in an array of in- 
creasing points ordered from left to right (xo, Xi, X2, x n , x n+ i) with x = x and 
x n+ i = 1 being the upper edge of the integration region. One constructs a rescaled 
array (x, x/x n , x/x 2 , x/xi, 1). We define = xjx^ and s n+ i = x < s n < s n _i < 
...Si < s = 1. We get 



At this point we introduce the linear interpolation 



(56) 



A(y) = 1 - V X _ l I A( Xi ) + V X _ l A(x t+1 ) (57) 
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and perform the integration on each subinterval with a change of variable y— > x/y 
and replace the integral J(x) with its discrete approximation Jn(x) to get 

J N (x) = A(x )^— f 1 ^-P(y)(y - Sl ) 

N 



1 - si J Sl y 

r dy 



+ E^WrVr ^P(y)(y-s i+1 ) 



i=l 

AT 



I 1 <7// 



^Afc) — — f" 1 ^P(y)(y - 
j = i — Sj y 



(58) 



with the condition j4.(xjv+i) = 0. Introducing the coefficients W(x,x) and x), 
the integral is cast in the form 



J N (x) = W(x, x)A(x) + Wfa, x)A(xi) 



i=i 



where 



W(X,X) = -J— ['^-(y-sJPiy), 
1 - Si 7 S i y 

^(Xi.x) = — *— /* ^ (y - s m ) P(y) 

rfy 



s i J si y 



(y-s,_!)P(y). 



We recall that 



f dx /<£)=/ ^ 

Jo (1 — #)+ JO 1 



/(y)-/(i) 



and that 



(1-x) 



1 ^ f ( \ — f 1 d v yf(y) - x f( x ) . -7 x , n x 

® /W = / + /W lo g(! - x ) 



x y y-x 



as can can be shown quite straightforwardly. 
We also introduce the expressions 



(59) 



(60) 
(61) 

(62) 



In (x) 
Jrn(x) 



1-S! 
1 



S 7 ; s< 



log(si) + log(l - Si) 



log ( + s i+1 log I 1 - * +1 



+1 



1 - Si 



1 ~~ s i 
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s i-l I ) 



1 &i 

Jnti(x) = logsi. 

1 - Si 



Using the linear interpolation formula (|57D we get the relation 

dyyA n (y) - xA n (x) . . . . 

= -log(l -z)A n (a;) + A n (x)/n (a;) 

* y y - a? 

jV 

+ ^ A n (xj) (Jrii(x) - Jnti(x)) (64) 
i=i 

which has been used for a fast and accurate numerical implementation of the 
recursion relations. 
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Figure 2: Fragmentation function for the up (top graph) and down quark (bottom 
graph) xD™(x, Q 2 ) at the lowest scale tinput) Qo = 10 GeV, and their QCD (or 
regular) and SQCD/QCD evolutions with Qf =10 3 GeV.The susy fragmentation 
scale is chosen to be 200 GeV. 
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Figure 3: Fragmentation functions into protons of the charm and strange quarks 
xDP'P(x,Q 2 ) at the lowest scale (input) Qo = 10 GeV, and its evolved QCD (regu- 
lar) and SQCD/QCD evolutions with Qf =10 3 GeV.The susy fragmentation scale is 
chosen to be 200 GeV. 




Figure 4: The gluon fragmentation function xD^(x, Q 2 ) at the lowest scale (input) 
Qo = 10 GeV, and its evolved QCD (regular) and SQCD/QCD evolutions with 
Qf =10 3 GeV.The susy fragmentation scale is chosen to be 200 GeV. 
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Figure 5: The fragmentation functions of squarks and gluino at the lowest scale 
(input) Qq = 10 GeV, with Q f = 10 3 GeV and susy scale 200 GeV. 
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Figure 6: The up quark fragmentation function xD^' p (x, Q 2 ) evolved with 
SQCD/QCD for varying final values of Qf. The susy fragmentation scale is cho- 



sen to be 200 GeV 
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Figure 7: The gluon fragmentation function xD^(x, Q 2 ) evolved regularly and ac- 
cording to SQCD/QCD for varying final values of Qf. The susy fragmentation scale 
is chosen to be 200 GeV. 
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Figure 8: The up quark fragmentation function xD^(x, Q 2 ) evolved regularly and 
according to SQCD/QCD for varying final values of the susy scale m 2 A and Qf = 10 s 
GeV. The susy fragmentation scale is chosen to be 200 GeV. 
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Figure 9: The quark and squark (up) fragmentation functions for a varying ni\ and 
a fixed large final scale Qf = 10 7 GeV. 
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Figure 10: The gluino fragmentation function xD p,p (x, Q 2 ) evolved in the region 
500 — Qf GeV for different values of Qf The susy fragmentation scale is chosen to be 
500 GeV. 
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